clear


% Choose lags, controls, etc.
lags = 4;
Tirf = 12;
trend = 0;          %=1 if detrend data (linear).
firstdiff = 0;      %=1 if APL is first-differenced.
uselagflows =0;     %=1 if first stage (Y/L) includes vars other than u or f.




% Read in Model IRFs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Choose calibration
%
baseline = 1;
if baseline ==1
    printinactionrate = 22.5;
    mat = readmatrix('/home/Output/Chains.xlsx','Sheet','Figure 8','Range','A2:J1003');
else
    % NOTE: When you run the C=0 case, you will see "Warning: Matrix is
    % close to singular". This happens for specifications that include
    % both lags of APL and lags of the JF rate. These are NOT the 
    % models we report in the paper. But they are included here
    % for completneness.
    printinactionrate = 0;
    mat = readmatrix('/home/Output/Chains.xlsx','Sheet','transition_OJS','Range','A2:J1003');
end
%
mo = mat(:,1);
Y_neg = mat(:,3);   % "neg" means impulse is -1%.
U_neg = mat(:,5);
F_neg = mat(:,6);
EU_neg= mat(:,7);
EE_neg= mat(:,8);
Hire_neg = mat(:,9) ./ (1-U_neg);
Newchain_neg = mat(:,10);
Lenchain_neg = mat(:,9) ./ Newchain_neg;

% Extract quarterly data
qt = (3:3:60)';
qts= length(qt);
diff = mo*ones(1,qts) - ones(length(mo),1)*qt';
[~,pos] = min(abs(diff));
pos = pos';
%
yss = Y_neg(1);
uss = U_neg(1);
fss = F_neg(1);
euss= EU_neg(1);
eess= EE_neg(1);
hiress = Hire_neg(1);
chainss = Lenchain_neg(1);
pss = yss / (1-uss);
%
y_mod = Y_neg([2;pos]);
u_mod = U_neg([2;pos]);
f_mod = F_neg([2;pos]);
eu_mod= EU_neg([2;pos]);
ee_mod= EE_neg([2;pos]);
hire_mod= Hire_neg([2;pos]);
chain_mod= Lenchain_neg([2;pos]);
%
logy = log(y_mod);
logp = log(y_mod) - log(1-u_mod);
logu = log(u_mod);
logf = log(f_mod);
logeu= log(eu_mod);
logee= log(ee_mod);
loghire = log(hire_mod);
logchain= log(chain_mod);




% IRF Matching
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Will match empirical IRF of APL. Read in:
pirf_data = readmatrix('/home/Output/IRFs.xlsx','Sheet','Fig 9','Range','B3:B14'); % (p,u,f) system.
%
% Simulation length
T    = 25000;
obs  = T - lags;
%
% Innovation std dev.
% These are chosen to match the empirical standard deviation of 
% HP filtered log APL = 0.01489 (over period 1994-2016).
sigma= .006208*baseline + .007158*(1-baseline);
%
% Initialize rho:
rho0 = .93*baseline + .91*(1-baseline);  
%
rng(502);
tfpeps = sigma*randn(T+1,1);    % innovations to aggregate TFP.
% 
% model vars and derivative
shock   = -.01;
logvar  = [logp, logu, logf, logeu, logee, loghire, logchain];  
ssvar   = [pss,  uss,  fss,  euss,  eess,  hiress,  chainss];
simvars = length(ssvar);
Dvar    = (logvar(2:end,:) - logvar(1:end-1,:)) / shock;  % "derivative" of outcome in time t w.r.t. innovation in t=0.
%
residtol = Tirf*(.0001^2);      % -> Avg squared difference bewteen model and data is .0001^2;
dresidtol= 1e-6;                % If percent decline in residual error < tolerance, then end. 1e-4
residlag= 1;
dresid= 1;
maxsim= 100;
resid= 10;
sim= 1;
%
format long g
while residtol < resid && dresidtol <dresid && sim <=maxsim 

    % Execute 3 times - one for initial rho, one for rho +.01 and one for rho -.01
    Drho = [0; .01; -.01];
    %
    pirf2_mod = zeros(Tirf, simvars -2, 3);   %add 1 to simvars b/c simvars omits quits (since quits=EE in model).       
    uirf2_mod = zeros(Tirf, simvars -2, 3);
    yirf2_mod = zeros(Tirf, simvars -2, 3);
    pirf3_mod = zeros(Tirf, simvars -3, 3);
    uirf3_mod = zeros(Tirf, simvars -3, 3);
    firf3_mod = zeros(Tirf, simvars -3, 3);
    yirf3_mod = zeros(Tirf, simvars -3, 3);
    %
    uirf2z_mod= zeros(Tirf, simvars -2, 3);
    yirf2z_mod= zeros(Tirf, simvars -2, 3);
    uirf3z_mod= zeros(Tirf, simvars -3, 3);
    firf3z_mod= zeros(Tirf, simvars -3, 3);
    yirf3z_mod= zeros(Tirf, simvars -3, 3);
    %
    for k = 1:3 
        rho = rho0 + Drho(k);
        if rho >=.999
            rho = .5*(rho0 + .999);
            Drho(k) = rho - rho0;
        end
        %
        % Draw innovations:
        logtfp_lag = -.5*(sigma^2);
        logtfp= zeros(T+1,1);
        for t = 1:T+1
            logtfp(t) = -.5*(sigma^2) + rho*logtfp_lag + tfpeps(t);     % E[e(eps)] = 1
            logtfp_lag= logtfp(t);
        end
        z = logtfp(2:end) - logtfp(1:end-1);
        zz= [zeros(qts-1,1);  z];   % initial history of zeros.
        %
        % Simulate paths:
        dlogvar = zeros(1,simvars);
        dlogvar_sim = zeros(T,simvars);
        for t = 1:T
            zt = sortrows( [(qts:-1:1)', zz(t:t+qts-1)], 1);   
            zt = zt(:,2);
            % sort in order of most recent to least recent.
            % eg, if t=3, then zt = [eps(3); eps(2); eps(1); 0->0];
            % eps(3) is then multiplied by Dvar(1), the contemporanous
            % effect of today's innovation; eps(2) is multiplied by 
            % Dvar(2), the contemporaneous effect of last period's
            % innovation; and so on. 
            %
            dlogvar = dlogvar + zt'*Dvar;
            dlogvar_sim(t,:) = dlogvar;
        end
        %
        % Run IRFs:
        var_sim = repmat(ssvar,T,1).*exp(dlogvar_sim);
        p = log(var_sim(:,1));
        u = log(var_sim(:,2));
        f = log(var_sim(:,3));
        eu= log(var_sim(:,4));
        ee= log(var_sim(:,5));
        hire= log(var_sim(:,6));
        chain = log(var_sim(:,7));
        MainData = [p,u,f,eu,ee,hire,chain]; 
        [~,vars] = size(MainData);
        %
        irfregs_runols;
        zirfregs_runols;    % Redo regs with true TFP innovations.
        irfregs_makeirfs;
        zirfregs_makeirfs;  % Redo IRFs based on true TFP innovations.
        %
        pirf2_mod(:,:,k) = reshape( irf2(:,1,:), Tirf, varseps2);    %vareps2 defined in makeirfs.m
        uirf2_mod(:,:,k) = reshape( irf2(:,2,:), Tirf, varseps2);
        yirf2_mod(:,:,k) = reshape( irf2(:,3,:), Tirf, varseps2);
        pirf3_mod(:,:,k) = reshape( irf3(:,1,:), Tirf, varseps3);
        uirf3_mod(:,:,k) = reshape( irf3(:,2,:), Tirf, varseps3);
        firf3_mod(:,:,k) = reshape( irf3(:,3,:), Tirf, varseps3);
        yirf3_mod(:,:,k) = reshape( irf3(:,4,:), Tirf, varseps3);
        %
        uirf2z_mod(:,:,k)= reshape( irf2_z(:,1,:), Tirf, varseps2);
        yirf2z_mod(:,:,k)= reshape( irf2_z(:,2,:), Tirf, varseps2);
        uirf3z_mod(:,:,k)= reshape( irf3_z(:,1,:), Tirf, varseps3);
        firf3z_mod(:,:,k)= reshape( irf3_z(:,2,:), Tirf, varseps3);
        yirf3z_mod(:,:,k)= reshape( irf3_z(:,3,:), Tirf, varseps3);
        %
        if k ==1
            ModTimeSeries= MainData;   % Save for use below.
        end
    end
    %
    % Compute output per worker residual 
    % Select BASELINE specification to 
    % use for model evaluation. 
    pirf_mod  = pirf2_mod(:,1,1);
    resid = (pirf_data - pirf_mod)'*(pirf_data - pirf_mod);
    disp( [sim, resid] );
    %pause;
    dresid= (residlag - resid) / residlag;
    residlag = resid;
    sim = sim+1;
    %
    % Update guess:
    pirf_mod_h = pirf2_mod(:,1,2);  % higher rho
    pirf_mod_l = pirf2_mod(:,1,3);  % lower rho 
    Jp = (pirf_mod_h - pirf_mod_l) ./ (Drho(2) - Drho(3));
    diffrho = (Jp'*Jp) \ (Jp'*(pirf_data - pirf_mod));
    newrho  = rho0 + (1/3)*diffrho;     % set step size.
    newrho  = max(.005, min(.995, newrho));
    disp( [rho0, newrho] );
    %pause;
    if residtol < resid && dresidtol <dresid
        rho0 = newrho;
    end
    w = warning('query','last');
    if isempty(w) ==0
        warnid = w.identifier;
        warning('off',warnid);
    end
end
if isempty(w)==0
    warning('on',warnid);
end
%
%
%
% Relative standard deviation of simulated series:
StdTimeSeries = std(ModTimeSeries);
RelStdTimeSeries = StdTimeSeries(2:end) / StdTimeSeries(1); 
titles  = ["Unemployment", "Job finding rate", "EU rate", "Quit rate", "Hiring rate", "Chain length"];
disp([titles; RelStdTimeSeries]);
%
% Print to Excel file. These are the results reported in Table G.
if baseline==1
    writetable(array2table(RelStdTimeSeries, 'VariableNames', titles), '/home/Output/IRFs.xlsx', 'Sheet', 'Table G', 'Range', 'B1:G2');
else
    writetable(array2table(RelStdTimeSeries, 'VariableNames', titles), '/home/Output/IRFs.xlsx', 'Sheet', 'Table G', 'Range', 'B4:G5');
end
%
%
% Save Model IRFs 
%
% As we did for the empirical IRFs, we present APL and U/L based on the 
% the first VAR system, {APL, U/L, JF rate}. The quit rate is 
% pulled from the second VAR system, {APL, U/L, quit rate}.
% Likewise, chain length is pulled from the final VAR 
% system, {APL, U/L, chain length}.
%
modirfs = [uirf2_mod(:,1,1), yirf2_mod(:,1,1),   yirf2_mod(:,3,1),  yirf2_mod(:,5,1)];
titles  = ["Unemployment", "Job finding rate",  "Quit rate",  "Chain length"];
% Print to Excel file. This is Figure 9. 
if baseline==1
    writetable(array2table(modirfs, 'VariableNames', titles), '/home/Output/IRFs.xlsx', 'Sheet', 'Figure 9', 'Range', 'B2:E14');
    writematrix([rho0;sigma], '/home/Output/IRFs.xlsx', 'Sheet', 'Figure 9', 'Range', 'C18:C19');
else
    writetable(array2table(modirfs, 'VariableNames', titles), '/home/Output/IRFs.xlsx', 'Sheet', 'Figure 9', 'Range', 'G2:J14');
    writematrix([rho0;sigma], '/home/Output/IRFs.xlsx', 'Sheet', 'Figure 9', 'Range', 'H18:H19');
end        
%
%
% Results based on true p innovations (Appendix H).
%
modirfs = [uirf2z_mod(:,1,1), yirf2z_mod(:,1,1),  yirf2z_mod(:,3,1), yirf2z_mod(:,5,1)];
% Print to Excel file. These are the results reported in Figure H2.
if baseline==1
    writetable(array2table(modirfs, 'VariableNames', titles), '/home/Output/IRFs.xlsx', 'Sheet', 'Figure H2', 'Range', 'B2:E14');
else
    writetable(array2table(modirfs, 'VariableNames', titles), '/home/Output/IRFs.xlsx', 'Sheet', 'Figure H2', 'Range', 'G2:J14');
end        
%
% Connection between regression residuals and true p innovations:
tfpshocks     = tfpeps(lags+2 : end);
coef_eps_res2 = zeros(1, simvars-2);
Rsq_eps_res2  = zeros(1, simvars-2);
for i = 1:simvars-2 
    %
    XX = [peps2(:,i), ones(T-lags,1)];
    coef = (XX'*XX)\(XX'*tfpshocks);
    coef_eps_res2(i) = coef(1);
    %
    resid = tfpshocks - XX*coef;
    Rsq_eps_res2(i) = 1 - ((resid'*resid) / ( (tfpshocks-mean(tfpshocks))'*(tfpshocks-mean(tfpshocks)) ));
end
%
% This allows for specifications where lags other than APL and U/L
% enter in the first stage. But if uselagflows = 0, then the 
% first stage is invariant to the outcomes examined in the
% second stage. Therefore, it is sufficient to just take
% the first element of the coefficients and R^2 vectors.
%
regmat = [coef_eps_res2(1), Rsq_eps_res2(1)];
titles  = ["Coefficient", "R squared"];
if baseline==1
    writetable(array2table(regmat, 'VariableNames', titles), '/home/Output/IRFs.xlsx', 'Sheet', 'Figure H2 Regressions', 'Range', 'B1:C2');
else
    writetable(array2table(regmat, 'VariableNames', titles), '/home/Output/IRFs.xlsx', 'Sheet', 'Figure H2 Regressions', 'Range', 'B4:C5');
end 














